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1.  Introduction 


To  accurately  describe  the  propagation  of  a  light  beam  through  a  medium,  one  must  account, 
roughly  speaking,  for  three  phenomena:  diffraction,  refraction,  and  absorption.  Diffraction 
refers  to  the  tendency  of  wave  phenomena  such  as  light  or  sound  to  deviate  from  rectilinear 
propagation  when  a  portion  of  the  wave  front  is  obstructed  in  some  way.  Diffractive  effects  are 
present  to  some  degree  in  every  real  optical  system,  simply  because  no  real  system  can  be 
infinite  in  extent;  every  real  system  must  necessarily  contain  an  effective  aperture  of  some  sort. 
Refraction  is  related  to  the  speed  with  which  a  wave  travels  through  a  medium.  It,  too,  can  result 
in  the  deviation  of  a  beam  from  straight-line  propagation,  as  when  the  refractive  properties  of  the 
medium  depend  on  position  in  such  a  way  that  different  portions  of  the  beam  experience 
different  amounts  of  refraction.  Absorption  is  the  process  responsible  for  the  change  in 
amplitude  of  the  wave  as  it  propagates  through  the  medium  and  energy  is  transferred  from  the 
wave  to  the  medium  or  from  the  medium  to  the  wave  (gain).  Both  refraction  and  absorption 
depend  on  the  amplitude  of  the  propagating  wave.  In  certain  media,  this  dependence  is  highly 
nonlinear,  while  in  other  materials,  such  as  “slow”  (fluence  dependent)  excited  state  absorbers 
(1 ),  it  is,  to  be  sure,  linear  in  the  wave  amplitude  but  reflects  the  entire  history  of  the  interaction 
between  the  medium  and  the  propagating  beam.  In  general,  calculating  the  effects  of 
propagation  of  a  high-intensity  laser  beam  through  all  but  the  simplest  media  presents  a  very 
difficult  problem. 

In  these  “simple”  media  (isotropic  media  exhibiting  refractive  and  absorptive  responses  that  are 
both  linear  and  time  independent),  beam  propagation  essentially  reduces  to  pure  diffraction. 
However,  even  a  pure  diffraction  problem  in  an  uncomplicated  geometry  can  generally  not  be 
solved  without  recourse  to  various  approximations,  and  these  impose  limits  on  the  resulting 
solutions.  Nevertheless,  solutions  do  exist. 

These  solutions  are  very  valuable.  The  special  case  of  beam  propagation  (diffraction)  in  an 
isotropic,  linear,  and  time-independent  medium  is  extremely  important  in  practice.  What  is 
more,  the  solution  of  the  linear  diffraction  problem  can  be  employed  to  attack  the  more  general 
problem  of  beam  propagation  in  a  nonlinear  medium  through  the  use  of  a  numerical  procedure 
that  separately  computes  the  effects  of  diffraction,  nonlinear  refraction,  and  nonlinear  absorption 
for  each  propagation  step.  This  is  the  “split  step”  Fourier  method  introduced  by  Fleck,  Morris, 
and  Feit  (2)  (sometimes  referred  to  simply  as  the  Beam  Propagation  Method). 

The  remainder  of  this  report  is  organized  as  follows.  Section  2  opens  with  an  examination  of  the 
general  nonlinear  optical  wave  equation  and  the  conditions  under  which  it  reduces  to  the  scalar 
Helmholtz  equation.  This  is  followed  by  a  careful  derivation  of  the  paraxial  wave  equation. 
Section  3  briefly  reviews  the  salient  results  of  classical  scalar  diffraction  theory.  In  sections  4 
and  5,  standard  integral  transform  methods  are  used  to  obtain  a  general  solution  to  the  paraxial 
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wave  equation  in  a  linear  medium  and  the  Helmholtz  equation  in  a  linear  medium,  respectively. 
These  results  underscore  the  fact  (perhaps  not  widely  appreciated)  that  the  exact  Rayleigh- 
Sommerfeld  integral  is  equivalent  to  the  scalar  Helmholtz  equation,  while  the  Rayleigh- 
Sommerfeld  integral  in  the  Fresnel  approximation  is  equivalent  to  the  paraxial  wave  equation. 
The  general  results  of  section  4  are  specialized  to  situations  with  cylindrical  symmetry  in  section 
6.  In  section  7,  these  cylindrically  symmetric  results  are  applied  to  a  practical  example — a 
Gaussian  beam  clipped  by  an  aperture. 


2.  Scalar  Electromagnetic  Wave  Equations  and  Their  Validity 


2.1  General  Wave  Equation  for  Nonlinear  Optical  Media 

The  derivation  of  an  electromagnetic  wave  equation  from  the  Maxwell  equations  is  presented  in 
virtually  any  textbook  on  optics  or  electromagnetic  theory.  The  most  general  fonn  (3)  of  the 
wave  equation  in  a  nonmagnetic  medium  containing  neither  free  charges  nor  free  currents  is,  in 
Gaussian  units, 


Vx  V  xE(j,z,?)  + 


1  ^ 

— t^E(x,z,0 
c  dt 


4 n  d2 
~^~Dt2 


P(x,z,t). 


(1) 


(In  anticipation  of  future  developments,  we  place  one  of  the  coordinates,  z,  on  a  separate  footing 
and  combine  the  remaining  two  coordinates  into  the  two-component  vector  x ,  which  designates 
positions  in  the  transverse  plane.)  Inequation  1,  P  (x,z,t)  is  the  total  material  polarization.  It  is 
usually  expressed  as  a  power  series  in  the  electric  field  E (x,z,t) ,  with  its  z'th  component  given 
by 

oo 

P .(x,z,t)=  ^dt'  x?](x,z,t-t')  Ej(x,z,t')  + 

-oo 

oo  oo 

J dt"  J dt' Z% (x, z,t-t',t-  t")Ej (x, z, t')Ek (x, z, t")  + 

-oo  -oo 

oo  oo  oo 

J  dt"'  J  dt"  J  dt'  Xifti(^’Z,t  - 1'  ,t  ~t",t  -  tm)E  j(x,z,t')Ek(x,z,t")E,(x,z,t'")  + 


In  this  expression,  the  quantities  x\"\ \{x,z,t,...)  are  components  of  the  nth-order  time-dependent 

electric  susceptibility  tensor,  and  the  repeated  indices  j,  k,  /,  ....  carry  an  implied  sum  (Einstein 
summation  convention).  The  convolutions  in  each  term  of  the  above  expansion  for  P,(f,z,t) 

become  simple  products  upon  Fourier  transformation. 
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We  consider  now  the  single-frequency  component  of  the  electric  field  oscillating  at  angular 
frequency  co\ 

E(x,  z,  t )  =  E(x,  z)e~uot  +  c.c.  (2) 

We  employ  the  well-known  vector  identity  for  the  curl  of  a  curl  and  the  constitutive  relation  in 
the  form 

4 nP  =  (n  -1  )E 


to  recast  equation  1  as 

viv-E)~V2E-C^EE  =  ^T(n2(x,z,(o,E;t)-n20)E .  (3) 

c  c 

Here,  n(x,z,a>,E;t )  is  the  Fourier  component  of  the  refractive  index  of  the  medium  oscillating 
at  frequency  co.  Strictly  speaking,  a  quantity  in  the  frequency  domain  is  a  function  of  co  only,  not 
of  t.  In  the  present  case,  however,  we  write  n  with  an  explicit  time  dependence  to  emphasize  the 
fact  that  each  Fourier  component  may  slowly  change  with  time  (here,  “slowly”  means  that  the 
time  scale  of  the  variation  is  much  longer  than  ah1)  as  the  properties  of  the  medium  evolve  as  a 
result  of  the  interaction  with  the  beam.  The  value  of  the  refractive  index  before  the  arrival  of  the 
beam  is  denoted  by  no. 

If  the  V (V  ■  E)  term  can  be  dropped,  equation  3  takes  the  form  of  the  inhomogeneous  Helmholtz 
equation,  which  can  be  solved  numerically  with  the  “split  step”  Fourier  scheme  alluded  to  in  the 
previous  section.  Roughly  speaking,  the  split  step  procedure  treats  effects  arising  from  the  right- 
hand  side  of  equation  3  separately  from  those  associated  with  the  left-hand  side.  One  solves  the 
homogeneous  Helmholtz  equation  over  a  short  longitudinal  step  and  then  accounts  separately  for 
the  effects  over  that  step  of  nonlinear  (and/or  time-dependent  linear)  refraction  and  absorption. 
The  homogenous  Helmholtz  equation,  together  with  the  simpler  equation  to  which  it  reduces  in 
the  case  of  paraxial  beams,  is  the  subject  of  the  remainder  of  this  report. 

We  postpone  to  a  future  report  a  full  analysis,  using  the  methods  of  Fax,  Fouisell,  and  McKnight 
(4),  of  the  implications  of  dropping  the  V(V  •  E)  term  in  equation  3.  For  the  present,  we  observe 
only  that,  since  free  charge  is  assumed  to  be  absent,  it  follows  from  the  Maxwell  equation  for  the 
electric  displacement  that  V  •  (n2E)  vanishes,  or  that 

V-E  =  -\E-V(lnn). 
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2.2  Homogeneous  Helmholtz  Equation 

In  a  linear  medium,  the  divergence  of  E  vanishes.  If,  in  addition,  the  properties  of  the  medium 
are  independent  of  time,  then  the  right-hand  side  of  equation  3  vanishes  as  well  and  equation  3 
reduces  to 

(v2  +  k2^E(x,z)  =  0  .  (4) 

Here,  k  =  n0a>/  c  ,  the  wave  number  in  the  medium.  In  an  absorptive  medium,  no  (and  thus  k ) 
has  a  positive  imaginary  part.  Equation  4  is  the  homogeneous  vector  Helmholtz  equation. 

We  designate  the  polarization  direction  of  the  field  with  the  unit  vector  e  and  write 
E(x,z )  =  eE{x,z ) .  If  e  is  constant,  equation  4  reduces  to  a  scalar  wave  equation, 

(V2  +k2)E(x,z)  =  0,  (5) 

the  scalar  Helmholtz  equation.  We  now  write  the  electric  field  as  an  envelope  function  times  a 
plane  wave  propagating  in  the  z-direction:  E(x,z)  =  A(x,z )  ek: .  In  terms  of  the  envelope 
function  A(x,z ) ,  equation  5  becomes 

(5 2  +  2 ikd  __  +  V 2  )A{x,  z)  =  0 ,  (6) 

where  V2  is  the  transverse  Laplacian,  given  in  Cartesian  coordinates  by  d 2  +  5 2 . 

2.3  Paraxial  Wave  Equation 

At  this  point,  we  proceed  with  the  slowly  vaiying  envelope  approximation  (also  known  as  the 
paraxial  approximation)  by  introducing  suitable  “non-dimensionalized”  quantities.  The  problem 
itself  sets  a  natural  length  scale  xo  for  the  transverse  coordinates,  e.g.,  the  size  of  an  aperture  or 
the  spot  radius  of  a  focused  laser  beam.  We  use  this  quantity  to  define  the  scaled  transverse 
coordinates:  X  =  x/x0.  The  transverse  length  scale  xo  and  the  wave  number  k  now  detennine 
z 0,  the  longitudinal  length  scale  or  “Rayleigh  range,”*  by  z0  =kx\,  and  the  corresponding  scaled 
longitudinal  coordinate  is  given  by  Z  =  z/zo.  We  denote  by  <?the  ratio  of  the  transverse  and 
longitudinal  length  scales:  e  =  xolzo={kxo)~  .  (An  alternate  procedure  for  non-dimensionalizing 
the  problem,  employed  by  some  workers  (5),  simply  expresses  all  lengths  as  multiples  of  the 
wavelength  X,  or  equivalently,  as  multiples  of  the  inverse  wave  number  k  '.  This  seemingly  very 
natural  choice  is  not  the  most  convenient  for  our  purposes.) 

In  terms  of  these  new  dimensionless  coordinates,  equation  6  becomes 

(s2d2z  +  2/3z+V2)z(l,Z)  =  0.  (7) 


If  one  chooses  x0  to  be  the  spot  radius  at  focus  of  a  Gaussian  beam,  then  z0  is  the  confocal  parameter  of  the  beam,  equal  to 
twice  the  Rayleigh  range  of  the  beam  as  conventionally  defined. 
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In  equation  7,  the  validity  of  the  slowly  varying  envelope  approximation  is  completely  manifest: 
as  long  as  the  transverse  length  scale  xo  is  much  greater  than  the  wavelength,  s«  1  and  the 
approximation  holds.  Equation  7  then  reduces  to 

(2idz  +V2±)a(X,Z)  =  0,  (8) 

the  paraxial  wave  equation.  In  terms  of  the  original  “dimensionful”  coordinates  (x,  z) ,  the 
paraxial  wave  equation  reads 

(likd2  +V2±)A(x,z)  =  0  .  (9) 


3.  A  Precis  of  Classical  Scalar  Diffraction  Theory 


3.1  History 

The  intuitive  picture  presented  in  1678  by  Christian  Huygens  (6)  forms  the  basis  for  the  modern 
theory  of  diffraction.  According  to  what  is  now  known  as  the  Huygens-Fresnel  principle,  every 
point  on  a  wave  front  may  be  thought  of  as  the  source  of  spherical  secondary  wavelets  that 
propagate  through  space  with  the  same  speed  and  frequency  as  the  primary  wave;  the  wave  front 
at  any  later  time  is  simply  the  envelope  of  the  secondary  wavelets.  In  his  famous  prize  paper 
presented  to  the  Paris  Academy  in  1818,  Augustin  Jean  Fresnel  combined  Huygens’s  envelope 
construction  with  Thomas  Young’s  ideas  about  interference.  The  Huygens-Fresnel  principle  was 
placed  on  a  firm  mathematical  basis  in  the  work  of  Gustav  Kirchoff  (7).  Kirchoff  s  theory, 
despite  its  inconsistent  assumptions  ( 8 ,  9 )  agrees  extremely  well  with  experiment  in  all  practical 
limits.  It  remained  only  for  Sommerfeld  (10)  and  Rayleigh  to  remove  the  inconsistencies  from 
the  Kirchoff  theory.  The  mathematical  expression  of  the  Huygens-Fresnel  principle  takes  the 
form  of  an  integral  which  (after  a  series  of  simplifying  approximations  that  impose  certain 
limitations  on  the  relative  size  of  the  diffracting  aperture,  the  distance  from  the  aperture  to  the 
plane  of  observation,  the  size  of  the  observation  region  on  the  plane,  and  the  wavelength  of  the 
radiation)  may  be  evaluated  analytically  in  a  variety  of  simple  aperture  geometries  (e.g., 
rectangle,  circle,  and  infinite  one-dimensional  slit).  Conventional  textbook  treatments  (11,  12, 

13,  14)  typically  include  several  of  these  approximate  calculations. 

3.2  Rayleigh-Sommerfeld  Diffraction  Integral 

Although  excellent  discussions  of  the  conventional  integral  fonnulation  of  diffraction  theory  are 
given  elsewhere  (we  find  presentation  of  (14)  particularly  lucid),  in  this  section  we  briefly 
review  the  standard  treatment  in  order  to  emphasize  several  points  that  surface  in  our  subsequent 
discussion.  We  consider  the  situation  at  an  observation  point  P  with  coordinates  (x,  z)  ,  as 

shown  in  figure  1 .  (Throughout  this  discussion,  we  use  the  two-component  vector  x  to  indicate 
positions  in  the  transverse  plane.)  The  vector  r  with  components  (x -x',z)  points  from  an 
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Figure  1.  Diffraction  by  an  aperture:  geometry  of  the  problem. 


infinitesimal  source  at  coordinates  (x',0)  to  the  observation  point  P;  the  location  of  the  source 
point  ranges  over  the  aperture.  The  resulting  Rayleigh- Sommerfeld  diffraction  formula  is  then 


2  n 


ikr 

ff  d2x  E(x, 0)— 

J  J 

plane 

\r  J 

COS' 


(M). 


(10) 


where  r  =  |r |  =  ^jz2  +  |jc  -  x '|  and  represents  the  distance  between  the  infinitesimal  source  and 

the  observation  point,  while  k  is  the  wave  number,  given  in  tenns  of  the  wavelength  2  in  the 
medium  by  k  =  2/zt/L  The  integration  in  equation  10  is  formally  over  the  entire  xy  plane,  but  the 
Kirchoff  boundary  condition,  which  sets  E(x,Q)  =  0  for  points  outside  the  aperture,  restricts  the 

integral  to  the  two-dimensional  aperture  itself.  For  the  case  in  which  the  aperture  is  illuminated 
by  a  normally  incident  plane  wave,  it  has  recently  been  shown  that  equation  10  may  be  rewritten 
in  terms  of  a  one-dimensional  parametric  integral  around  the  perimeter  of  the  aperture  (15). 

The  direction  cosine  in  equation  10,  which  has  the  value  z/r,  is  known  as  the  obliquity  factor  or 
inclination  factor.  Fresnel  himself  recognized  the  need  for  a  quantity  of  this  kind  (12),  without 
which  the  secondary  wavelets  envisioned  by  Huygens  would  also  generate  a  reverse  wave 
propagating  back  toward  the  source  illuminating  the  aperture;  no  such  backward  propagating 
wave  is  observed  experimentally.  In  the  Fresnel-Kirchoff  diffraction  formula,  the  obliquity 
factor  is  }(cos(r,z)  + 1).  Interestingly  enough,  equation  10,  including  the  proper  form  of  the 

obliquity  factor,  can  be  derived  entirely  from  Fourier  transform  theory,  with  no  reference 
whatsoever  to  the  Huygens-Fresnel  principle  (5). 

Equation  10  treats  electromagnetic  radiation  as  a  scalar  phenomenon,  completely  ignoring  the 
fact  that  the  various  components  of  the  electric  and  magnetic  field  vectors  are  coupled  by 
Maxwell’s  equations  and  thus  cannot  be  treated  independently.  Within  the  limitations  of  a  scalar 
theory,  however,  equation  10  is  considered  to  be  exact  throughout  the  entire  diffraction  space. 
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Fortunately,  the  limitations  that  a  scalar  theory  places  on  the  macroscopic  length  scale  are  not 
particularly  severe;  experiments  involving  the  diffraction  of  microwaves  from  apertures  confirm 
that  the  scalar  theory  is  valid,  provided  the  aperture  size  is  much  greater  than  the  wavelength  of 
the  radiation  (16),  while  a  detailed  theoretical  analysis  involving  focused  Gaussian  laser  beams 
shows  that  the  scalar  treatment  suffices  so  long  as  the  focal  spot  size  is  at  least  approximately  as 
large  as  the  wavelength  (17). 

3.3  Preliminary  Approximations 

At  distances  z  from  the  aperture  on  the  order  of  the  wavelength,  the  exact  Rayleigh-Sommerfeld 
diffraction  fonnula,  equation  10,  can  be  evaluated  by  the  method  of  stationary  phase,  which  is 
described  in  (11).  For  z  values  of  at  least  several  wavelengths,  one  can  take  (l  /  r-  ik ) «  -ik  ,  in 

which  case  equation  10  becomes 

E(x,z)  =  —  ff  d2x'  E(x\ 0)ear4-  (H) 

A  r 

plane 

The  highly  oscillatory  character  of  the  integrand  makes  the  direct  numerical  evaluation  of 
equation  1 1  problematic;  however,  a  recent  report  suggests  that  present-day  desktop  computers 
running  commercially  available  applications  software  may  now  be  equal  to  the  task,  at  least  for 
configurations  for  which  the  Fresnel  number  is  not  too  large  (18). 

3.4  Fresnel  Approximation 

The  standard  method  of  attacking  the  integral  in  equation  1 1  relies  on  the  expansion 


'=!  +  - 


x  ■ x 


2  z2 


2z2 


|2  f 

■  +  o  ■ 


X -x 


v  z  7 


One  approximates  the  factor  of  r  in  the  denominator  of  the  integrand  with  only  the  zeroth-order 


tenn,  while  for  the  factor  of  r  in  the  phase  eh ,  the  three  quadratic  terms  are  retained  as  well. 

Jkz  r  ik 


E(x,z)  =  —~ —  ff  d2x  E(x',0)  exp 
A  z 


plane 


2  z 


\x  -x 


(12) 


This  is  the  Fresnel  approximation,  whose  validity  has  been  the  subject  of  a  number  of 
investigations  over  the  years  (19,  20). 

3.5  Fraunhofer  Limit 

If  now  the  dimensions  of  the  aperture  are  such  that  for  every  point  (x',0)  in  the  aperture, 

k  1x1" 
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2 

then  exp[4£  x'  /  z\  «  1  and  equation  12  reduces  to 


X  z 

the  Fraunhofer  approximation. 


•  i&Z 

E(x,z)  = - exp[|&  |%|2/z]  JJ  d2xf  £(x',0)exp 


plane 


z 


4.  Solution  of  the  Paraxial  Wave  Equation 


One  can  easily  verify  by  direct  substitution  that  equation  12,  the  diffraction  integral  in  the 
Fresnel  approximation,  solves  the  paraxial  wave  equation,  equation  9,  with  E(x,z)  =  A(x,z )  elkz  , 
thus  demonstrating  the  equivalence  of  the  Fresnel  and  paraxial  approximations.  It  is  nevertheless 
illuminating  to  proceed  with  the  direct  solution  of  equation  7  by  standard  methods,  imposing  the 
“paraxial  limit”  s  — »  0  only  when  necessary. 


4,1  Helmholtz  Equation  as  a  Singular  Perturbation  of  the  Paraxial  Wave  Equation 

We  express  A(X,Z)  in  terms  of  its  Fourier  transfonn  in  the  transverse  coordinates. 

A(X,Z)  =  (2k)  '  |J  d~K  A{k,Z )  eig'* 

plane 


The  Fourier- transformed  field  amplitude  A{k,Z)  satisfies 

(s2d2z+2idz-K2)A{K,Z)  =  0,  (13) 

where  k  =  |j?|  .  Setting  A(k,Z)  =  (//(/?)  e,aZ  ,  we  find  that  equation  13  is  satisfied,  provided  that 

s2 a2  +  2a  +  k2  =  0  .  (14) 


There  are  two  roots, 


«+ 


-1  +  Vl -s2k2 


-\k2  +o 


and 


+  \k^  +o(e 

Returning  to  the  physical  field, 

equation  2, 


—  1  —  VI  -£  K  -2 

a  - - ; - =  —2s 


of  which  the  second,  a_ ,  is  singular  in  the  limit  s  — >  0 . 
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E(X,Z,t)  =  e(2/r)  1  fjd2/c  i//(/c)  e,K'x  exp[z'(£  2  +  a±)Z  -ia>t]  +  c.c. 

plane 

=  e(2 ny1  JJ d2tc  y/(K)  e,K'x  exp[±z  -  icot]  +  c.c. 

plane 

we  see  that  the  root  a_  corresponds  to  a  solution  that  propagates  in  the  negative  z-direction. 

Nowhere  in  the  region  Z  >  0 ,  either  at  infinity  or  elsewhere,  is  there  a  source  from  which  such  a 
backward  propagating  wave  could  originate,  so  we  are  justified  in  discarding  a_  as  unphysical. 

Thus,  whereas  in  the  integral  formulation  of  diffraction  theory,  one  introduces  the  obliquity 
factor  in  order  to  suppress  troublesome,  backward  propagating  Huygens  wavelets,  here  one  must 
exorcise  the  demon  with  an  appeal  to  the  boundary  condition  at  Z  — »  +oo  ! 

In  the  language  of  differential  equations,  one  would  say  that  the  term  of  order  A  in  equation  7  is 
a  singular  perturbation.  Such  a  perturbation  can  completely  alter  the  character  of  the  solution.  In 
the  present  case,  “turning  on”  the  perturbation  introduces  a  second  root  in  equation  14  where 
there  was  only  one  before  and  so  in  some  sense,  it  is  not  unnatural  that  this  new  root  should  be 
singular.  It  corresponds  to  a  rapidly  oscillating  solution  of  the  differential  equation,  and  this 
solution  can  be  very  troublesome  indeed. 

4.2  Convolution  Kernel  for  the  Paraxial  Wave  Equation 

To  obtain  a  solution  to  the  paraxial  wave  equation,  we  discard  the  singular  root  and  approximate 
the  remaining,  nonsingular  root  by  the  zeroth-order  term  in  its  expansion  in  powers  of  s  : 

A(X,Z)  =  (2 ttY1  JJ d2K  \j/{K)e~kKlz  eiiiX  (15) 

plane 

The  function  (//(/?)  is  determined  by  the  boundary  condition  on  A  at  Z  =  0  . 


y/{ic)  =  A(k, 0)  =  (2 ttY1  JJ  d2X'  A{X\ 0)  e~iiix'  (16) 

plane 


Substituting  equation  16  into  15  and  reversing  the  order  of  integration,  we  obtain 

d2x 


A(X,Z)=  JJ  d2X'  A(X', 0)  JJ 


plane 


plane 


[2zr  ]2 


eXK'z  eiH  (X-X') 


(17) 


Recalling  that  k2  =  k 2  +  k2  ,  we  see  that  the  integral  over  the  /v  plane  in  equation  17  factors  into 
a  product  of  two  identical  one-dimensional  integrals  of  the  form 


Gx(x-x\  z)  = 


e-N2z  eiq(x-xf) 


i(X-X')2 


^jlniZ 


(18) 
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(19) 


Equation  17  thus  takes  the  form  of  a  two-dimensional  convolution, 

A(X,Z)  =  JJ  d2X'  A(X',  0)  G(X  -  X',  Z )  , 

plane 


in  which  the  kernel, 


G(X-X',  Z)  = 


X-X'\ 


2 


2  7tiZ 


=  Gi(x-x\  z)Gl(r-r,  z), 


(20) 


is  the  square  of  the  one-dimensional  convolution  kernel  in  equation  18.  Expressed  in  the 
“dimensionful”  coordinates  ( x ,  z ) ,  the  two-dimensional  convolution  kernel  in  equation  20  is 


G(x  -x',  z )  = 


•  2 

«0 

Zz 


exp 


ik  i_  _,i2 

— \x  -x 
2z 


(21) 


Recalling  that  E(x,z)  =  A(x,z)  e'kz ,  we  see  that  equation  19  is  identical  to  equation  12,  the 
Fresnel  approximation  of  the  Rayleigh-Sommerfeld  diffraction  integral. 

4.3  A  Remark  on  Convergence 

In  order  to  ensure  convergence  of  the  integral  in  equation  18,  Z  =  z/  kx\  must  have  a  negative 

imaginary  part.  Since  z  and  xo  are  both  real,  the  requirement  that  Im[Z]  <  0  is  equivalent  to  the 
condition  Im[k]  =  a/2  >  0,  in  which  a  is  the  absorption  coefficient  of  the  medium;  mathematical 
convergence  is  thus  assured  in  the  case  of  propagation  through  an  absorptive  medium.  Of 
course,  the  integral  in  equation  1 8  also  converges  for  propagation  through  a  lossless  medium.  In 
this  case,  we  observe  that  a  perfectly  monochromatic  field  is  a  mathematical  idealization.  An 
actual,  physical  field  must  have  been  turned  on  at  some  time  in  the  past,  which  one  may  account 
for  by  adding  a  small  positive  imaginary  part  to  the  frequency,  so  that  e~,at,  and  thus  the  physical 
field,  vanishes  as  t  — >  -qo  .  This  is  equivalent  to  a  positive  absorption  coefficient. 


4.4  Numerical  Computation  of  the  Fresnel  Convolution  Integral 

In  practice,  the  integral  on  the  right-hand  side  of  equation  19  will  generally  have  to  be  evaluated 
numerically.  This  is  most  efficiently  accomplished  through  the  use  of  Fourier  transfonn  methods 
using  the  well  known  fast  Fourier  transform  (FFT)  algorithm.  Given  a  complex-valued  function 
/ ( X )  on  the  plane  9? 2 ,  we  symbolically  denote  the  two-dimensional  Fourier  transform  of 
/ ( X )  ,  which  is  a  function  of  k  ,  by 

Firm  =  f(K)  =  [2 k\x  JJ d2X  f(X)e-*  * 

plane 

and  the  inverse  transform  of  / (/?) ,  a  function  of  X  ,  by 
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F-'[/](X)  =  f(X)  =  \ln]1  JJ  d2K  f{ic )  e**  . 

plane 

Using  these  relations  and  the  well-known  integral  representation  of  the  Dirac  delta  function  in 
terms  of  complex  exponentials,  one  can  write  equation  19  symbolically  as 

A(X,Z)  =  2 n  F-‘[FM](£,  0)  F[G](/f,  Z)].  (22) 

Equation  22  is  a  prescription  for  calculating  A(X,Z) :  compute  the  Fourier  transform  of 
A(X,  0) ,  the  field  at  the  boundary,  Z  —  0;  multiply  this  by  the  Fourier  transform  of  the 
convolution  kernel,  equation  20;  compute  the  inverse  Fourier  transfonn  of  the  resulting  product. 
Now,  the  Fourier  transform  of  the  convolution  kernel  in  equation  20  is  given  by 

F [G](ic,  Z)  =  [27r]-le-^2/z, 


so  equation  22  is 


A(X,Z)  =  F  1 


F [A]{k,  0)  e41"1" 


This  is  precisely  equation  15,  in  which  y/{^)  ■>  given  by  equation  16,  is  simply  F [A~\{k,  0) ,  the 
Fourier  transform  of  the  boundary  condition  at  Z  =  0.  Thus,  we  see  that  same  procedure  that  we 
employed  to  derive  the  solution  in  equation  19  and  the  corresponding  Fresnel  kernel  in  equation 
20  also  provides  the  most  efficient  means  of  evaluating  the  integral  expression  in  which  the 
solution  is  given. 


5.  Solution  of  the  Helmholtz  Equation 


In  this  section,  we  derive  the  Rayleigh-Sommerfeld  result  by  employing  the  integral  transform 
methods  to  solve  equation  5,  the  scalar  Helmholtz  equation.  The  calculation  proceeds  exactly  as 
in  the  previous  section.  We  take 

E(x,z )  =  (2#)  1  JJ  d2K  y/(i<)  e'K  X  e'az ,  (23) 

plane 

where  is  detennined  by  the  boundary  condition  on  E  at  z  =  0: 

y/(K)  =  (2^  JJ  d2x  E{x, 0)  e**' .  (24) 

plane 

The  Helmholtz  equation  is  satisfied,  provided  that 

a  =  ±y/k2  -rc2  .  (25) 


11 


The  negative  root  corresponds  to  a  backward  propagating  solution  and,  as  before,  we  discard  it 
as  unphysical.  Substituting  equations  24  and  25  into  23  and  reversing  the  order  of  integration, 
we  obtain 

E{x,z)  =  (IttY2  JJ  d2x'  E(x', 0)  JJ  d2KeiH%-^  e ,  (26) 

plane  plane 


or,  expressed  in  the  dimensionless  variables  introduced  in  section  2.3, 


E(X,Z)  =  (2tvY2  JJ  d2X'E(X',  0>JJ  d~Ke 


2  iii(X-X')  izR-e2*1  Is 2 


plane 


plane 


Utilizing  the  result  of  Lalor  (27), 

7  2  //c-x  7zV k2  -K 


ii“ 

plane 


Ke  e  =  -2n — 

dz 


z  +  \x\ 


we  perform  the  integration  over  ron  the  right-hand  side  of  equation  26  and  so  obtain 


Her 

E(x,z)  =  ( 2k)  1  JJ  d2x'  E(x', 0)- — 


plane 


—  ~ik 
\r  J 


which  is  the  Rayleigh-Sommerfeld  formula,  equation  10,  with  obliquity  factor  z/r. 


6.  Axial  Symmetry 


We  now  turn  our  attention  to  problems  exhibiting  cylindrical  symmetry  about  the  z-axis.  As 
before,  xo  denotes  the  transverse  length  scale  characterizing  the  problem.  We  employ  this  length 
scale  to  define  the  dimensionless  radial  coordinate  R  =  r/x o  and  we  use  the  resulting  polar 
coordinates  ( R ,  6)  to  designate  positions  in  the  transverse  plane. 


6.1  Cylindrically  Symmetric  Restriction  of  the  General  Fresnel  Kernel 

We  impose  the  restriction  of  cylindrical  symmetry  on  equation  19,  the  general  solution  of  the 
paraxial  wave  equation,  and  on  the  associated  convolution  kernel  in  equation  20.  Cylindrical 
symmetry  requires  that  the  field  amplitude  A(R,Z)  be  independent  of  the  angular  coordinate,  so 
without  loss  of  generality,  we  set  6=0.  We  use  the  Law  of  Cosines  to  write 


X  -  X' 


=  R  -2RR'  cos  6'  +  R'  and  so  express  the  convolution  kernel  in  equation  20  as 


G{R,9,R'  ,6'  ,Z)\  =— !— expl 

IniZ 


2  Z 


(f?2 -2RR'cos6'  +  R'2) 
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Inserting  the  above  expression  for  the  convolution  kernel  into  equation  19  with 

d2X'  =  R'dR'dO' ,  we  employ  the  well-known  integral  representation  of  the  Bessel  function, 


2tt 

—  f  d&t 

2x{ 


-iu  cos  O' 


=  J0{u). 


to  perfonn  the  integration  over  9'  and  obtain  the  following  expression  for  the  field  amplitude. 

RR 


A(R,Z)  =  -—  exp 


—  R2 

2  Z 


jdR’R'AiR’,  0)e2Z  J{ 


r'2  (  DD'\ 

1 

0 


V  Z  j 


(27) 


We  write  equation  27  in  the  following  form: 

oo 

A(R,Z)  =  J dR'R'A(R',0)  G(R,R',Z) , 


(28) 


in  which  the  factor  of  R'  appears  explicitly  in  the  integrand.  The  associated  kernel  is  then 


G(R,R',Z)  =  —  exp 
iZ 


i(R~  +R'Z ) 
2Z 


Jr 


f  rrx 


(29) 


V  a  J 


the  analog  of  equation  20  for  situations  involving  cylindrical  symmetry.  (By  including  the  factor 
of  R'  in  the  integration  measure  instead  of  the  kernel,  we  obtain  an  expression  for  the  kernel  that 
is  symmetric  in  the  variables  R  and  R' .) 

6.2  Solution  of  the  Paraxial  Wave  Equation  Via  Hankel  Transform 

It  is  also  possible  to  obtain  the  cylindrically  symmetry  kernel,  equation  29,  in  a  somewhat  more 
methodical  fashion,  namely,  by  direct  solution  of  the  paraxial  wave  equation  in  two  transverse 
dimensions  under  the  assumption  of  cylindrical  symmetry.  In  this  case,  the  transverse  Laplacian 
operator  is  given  by  d\+ \d R  and  the  paraxial  wave  equation  8  takes  the  form 

{2idz+d\+idR)A{R,Z)  =  0. 

The  cylindrical  symmetry  of  the  problem  suggests  that  we  express  A(R,Z)  in  tenns  of  its  Hankel 
transform, 

oo 

A(R,Z)  =  jd/c  kA(k,Z)  J0{kR)  .  (30) 

o 

The  field  amplitude  in  the  transform  domain,  A(k,Z ) ,  satisfies  the  ordinary  differential  equation 
(; idz  -\k2)  A(k,Z)  =  0.  The  solution  is: 

~  ~  --k2Z 

A(fc,Z)=A(fc,0)e  2  ,  (31) 
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where  A(/c,  0)  is  given  by  the  inverse  Hankel  transform  of  the  boundary  condition  at  Z  =  0. 

oo 

A(k,0)  =  \dR'R'  A(R'$)J0(kR')  (32) 

o 

If  we  substitute  equation  32  into  31,  then  substitute  the  result  in  equation  30  and  reverse  the 
order  of  integration,  we  obtain 

co  co  --k2Z 

A(R,Z)  =  J dR'R'  A(R',0)  \dx  k  J 0(kR')  J 0(hR)e~ 2"  . 

0  0 

Comparing  this  last  expression  with  equation  28,  we  recognize  that  the  k  integral  is  none  other 
than  G(R,R',Z ) ,  the  kernel  that  we  seek: 

' RR * 

CZ~/ 


G(R,R',Z)  =  J  die  k  J0(kR')  J0(hR)e 


-  kZ 
2 


= - e 

Z 


l  — (R'-+R2) 

1  2  Z  T 

•J  n 


The  integral  is  given  in  Watson  (22)  and  reproduces  our  previous  result,  equation  29. 


6.3  Cylindrically  Symmetric  Fraunhofer  Kernel 

In  the  case  of  circular  symmetry,  the  Fraunhofer  condition  expressed  in  dimensionless  variables 
reads  R'2  /  (2Z)  «  1 .  In  the  Fraunhofer  limit,  the  factor  cxp[f /C2  / Z]  in  the  integrand  of 

equation  27  may  be  replaced  by  unity,  which  gives 


A(R,Z) 


exp 


2Z 


-R- 


^dRR'  A(R',0)J0 
0 


'RR* 


(33) 


Thus,  the  Fraunhofer  kernel  associated  with  equation  28  is  given  by 


G(R,R’,Z) 


1 

iZ 


exp 


iR2 
2  Z 


J n 


'RR* 


7.  Example:  Clipped  Gaussian  Beam 


We  conclude  with  an  example  that  is  of  considerable  practical  importance:  the  diffraction  of  a 
principal-mode  Gaussian  beam  by  a  coaxial  circular  aperture  that,  for  simplicity,  will  be  assumed 
to  lie  in  the  focal  plane  of  the  beam.  This  problem  is  characterized  by  two  radial  length  scales: 
jco,  the  radius  of  the  aperture,  and  wq,  the  half-width  of  the  beam  at  He ~  of  its  on-axis  intensity. 
Either  quantity  could  be  used  to  define  a  dimensionless  radial  coordinate;  we  choose  the  former, 
setting  R  =  r/x o.  We  also  define  Wo  =  m/xo,  the  dimensionless  quantity  characterizing  the  extent 
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of  the  beam.  The  clipped  Gaussian  beam  is  then  described  by  the  following  Kirchoff  boundary 
condition  at  Z=  0. 


\Ane~R/wi  R<  1 

A(R,0)=  0 

0,  R>  1 


(34) 


7.1  Fresnel  Diffraction  of  a  Clipped  Gaussian  Beam 

With  the  clipped  Gaussian  boundary  condition  in  equation  34,  equation  27  becomes 

1  f  RR ^ 


Z(A,Z)  =  -l^-eiR2/2Z\dR'R' 


-nR'j, 


\  £  J 


(35) 


where  77  =  W02  -i  /(2Z)  .  After  repeated  partial  integration  and  application  of  the  well-known 
recurrence  relation 


d_ 

dx 


[xnJn{x)\=x"J„^{x), 


equation  35  becomes 


A(R,Z)  =  —lA° 


2  r/Z 


exp 


iR 2 
2  Z 


V 


z 


'2/yZY 

A 


2 


(36) 


Associated  with  the  field  amplitude  of  equation  36  is  the  intensity  profile 

|2 


/(0,Z) 


=  e’-l 


1-2 


Z| 

n= 1 


'2^zy J  ( R 3 

R 


J 


(37) 


7.2  Fraunhofer  Diffraction  of  a  Clipped  Gaussian  Beam 

Substituting  equation  34,  the  boundary  condition  for  a  clipped  Gaussian  beam,  into  the 
Fraunhofer  result,  equation  33,  one  still  obtains  equation  35,  except  that  77  =  W02 ,  instead  of 

W()  2  -  i  /(2Z)  .  We  thus  obtain  the  field  amplitude  of  a  clipped  Gaussian  beam  in  the  Fraunhofer 
limit  by  simply  setting  r]  =  W()  2  in  our  previous  result,  equation  36. 

7.3  Limiting  Case:  Top  Hat  Beam  (Uniformly  Illuminated  Aperture) 

A  beam  that  overfills  the  aperture  (W0  »  1 )  provides  illumination  that  is  effectively  uniform.  In 
this  case,  77  «  -i /(2Z)  and  equation  36,  the  solution  in  the  Fresnel  regime,  reduces  to 

—  («2+l)  “  /  N.,,  ,  v 

A{R,Z)  =  \e*  lW*'z)' 
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As  we  saw  in  section  7.2,  the  Fraunhofer  limit  is  obtained  by  setting  77  =  W02  in  equation  36.  In 
the  present  case,  we  have  W0  »  1 ,  so  to  obtain  the  Fraunhofer  limit,  we  need  only  set  77  =  0  in 
equation  36.  Doing  so  yields 


A(R,Z) 


- R2 

-iAn  e2Z 


J,{R!Z) 

R 


The  associated  intensity  profile  in  the  far  field  is  then  given  by 


I{R,Z )  _  f  2  J^R/Zh2 

/(0,z)_t  R/z  , 


(38) 


the  77  — >  0  limit  of  equation  37.  Noting  that  R/  Z  =  kx0r  /  z  =  kx0  tana  «  kx0  sin  a  ,  where  a  is 

the  angle  that  the  observation  point  makes  with  the  z-axis  as  measured  from  the  center  of  the 
aperture,  we  recognize  that  equation  38  is  in  fact  the  celebrated  result  of  Airy  (23). 


8.  Summary 


From  the  Maxwell  equations  in  a  nonmagnetic,  isotropic  medium,  one  may  derive  a  general 
nonlinear  vector  wave  equation.  If  V(V  •  E)  =  0 ,  this  reduces  to  the  inhomogeneous  vector 
Helmholtz  equation,  which  is  generally  written  with  the  nonlinearity  contained  in  the  source  tenn 
on  the  right-hand  side.  If  the  direction  of  the  electric  field  polarization  is  constant,  as  is  usually 
the  case,  one  solves  the  corresponding  scalar  equation.  In  a  geometry  in  which  the  transverse 
length  scale  is  much  larger  than  the  wavelength  of  the  radiation,  the  Helmholtz  equation  is  well 
approximated  by  the  paraxial  wave  equation,  which  is  only  first  order  in  the  longitudinal  variable 
and  thus  much  easier  to  solve.  The  Helmholtz  equation  may  be  regarded  as  a  singular 
perturbation  of  the  paraxial  wave  equation,  and  some  of  the  difficulties  arising  in  the  solution  of 
the  former  partial  differential  equation  are  related  to  this  fact. 

The  Huygens-Fresnel  principle  provides  a  very  intuitive  picture  of  the  diffraction  of  light  from 
an  aperture.  The  mathematical  expression  of  the  Huygens-Fresnel  principle  takes  the  form  of  an 
integral  over  the  aperture:  the  Rayleigh-Sommerfeld  diffraction  integral.  The  oscillatory  nature 
of  the  integrand  makes  exact  evaluation  of  the  Rayleigh-Sommerfeld  integral  problematic,  so 
one  frequently  employs  the  simplifying  Fresnel  approximation  in  evaluating  the  integral.  We 
have  demonstrated  how  these  same  results  may  be  obtained  without  appeal  to  the  Huygens- 
Fresnel  principle  by  using  standard  methods  to  solve  the  appropriate  partial  differential 
equations:  the  exact  Rayleigh-Sommerfeld  integral  may  be  obtained  from  the  scalar  Helmholtz 
equation  in  a  linear  medium,  while  the  Rayleigh-Sommerfeld  integral  in  the  Fresnel 
approximation  results  from  solution  of  the  paraxial  wave  equation  in  a  linear  medium. 
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One  can  employ  the  solutions  of  these  linear  diffraction  problems  to  attack  the  more  general 
problem  of  beam  propagation  in  a  nonlinear  medium  through  the  use  of  the  numerical  procedure 
introduced  by  Feit  and  Fleck  in  the  late  1970s.  This  procedure,  known  as  the  “split  step”  Fourier 
method  or  simply  as  the  Beam  Propagation  Method,  separately  computes  the  effects  of 
diffraction,  nonlinear  refraction,  and  nonlinear  absorption  for  each  propagation  step. 
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